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We report parabolic bursting in a globally coupled network of mixed population of oscillatory 
and excitable Josephson junctions. The resistive-capacitive shunted junction (RCSJ) model of the 
superconducitng device is used for this study. We focus on the parameter regime of the junction 
where its dynamics is governed by the saddle-node on invariant circle (SNIC) bifurcation. In this 
SNIC regime, the bursting appears in a broad paramater space of the ensemble of mixed junctions. 
For a coupling value above a threshold, the network splits into two synchronized clusters when a 
reductionism approach is applied to reproduce the bursting behavior of the large network. The ex¬ 
citable junctions effectively induces a slow dynamics in the network to generate bursting. This 
bursting is a generic property of a globally coupled network with a mixed population of dynamical 
nodes where each node posseses the SNIC property. 

PACS numbers: 05.45.Xt, 05.45.Gg, 85.25.Cp, 87.19.1m 


I. INTRODUCTION 

The superconducting Josephson junction shows self¬ 
oscillation m for an applied constant current above 
a critical value. It is usually modeled as a resistive- 
capacitive-shunted junction (RCSJ) which has its me¬ 
chanical analog in a damped pendulum with a constant 
torque. A RCL-shunted junction (resistive-capacitive- 
inductively-shunted junction) model [4. 8] was also used 
to include an inductive loading effect in an array of junc¬ 
tions where more complex dynamics including chaos was 
seen. Interestingly, the superconducting device shows 
some typical spiking and bursting behaviors mm most 
commonly seen in a Type I excitability neuron m ■ The 
bursting dynamics was also found prominent mmm in 
a periodically forced junction. This is due to the intrinsic 
SNIC characteristic of the junction in a selected param¬ 
eter space Enana, which typically governs a class of 
bursting dynamics in Type I excitability neurons. 

Spiking is a repeatative firing state and bursting is a 
state of recurrent switching between a firing state or os¬ 
cillatory state and a resting state. The minimal condition 
for bursting in a system necessitates the presence of an 
intrinsic slow-fast dynamics mmm- As example, in 
biological neurons, the simplest ionic processes involved 
in spiking are due to the flow of Na+ and K + ions across 
the cell membrane, while the bursting may be observed 
when the fast spiking (FS) is controlled by a slow pro¬ 
cess like Ca ++ -gated K + ion movement across the mem¬ 
brane. The slow dynamics controls the firing or start of 
the oscillation and intermittently stops it when the tra¬ 
jectory of the dynamics moves slowly towards a steady 
state. Alternatively, an excitable system when coupled 
to an oscillatory system, was found m to induce a slow 
dynamics and thereby originates a type of chaotic burst¬ 
ing. 

On a different context, a mixed population of glob¬ 
ally coupled inactive or excitable and active or oscilla¬ 


tory units was investigated earlier 18U22] in search of 
synchrony and global oscillation. Such a global oscilla¬ 
tion is practically important, particularly, in the context 
of a desired synchrony of the pacemaker cells [23] [24]. It 
is also important to know, in the event of a growing cell 
death, how robust are the pacemaker cells in the heart or 
the suprachiasmatic cells in the brain to sustain a glob¬ 
ally synchronized oscillation? In the dynamical sense, a 
death of a cell is considered as a passive or an excitable 
state. In the situation of a progressive cell death, in other 
words, increasing number of passive oscillators, a popula¬ 
tion of globally coupled oscillators showed a type of aging 
transition [IS]. Such aging transition or death state is not 
the focus of this current study. We emphasize rather on 
the synchronized state (1:1 or higher phase-locking) of 
global oscillation of the mixed population as shown ear¬ 
lier [18 ] 19] where the type of oscillatory dynamics was 
not given appropriate attention. 

In this backdrop, we consider the superconducting 
RCSJ model to construct a globally coupled network of 
mixed population of excitable and oscillatory junctions 
and, particularly, focus on its collective coherent dynam¬ 
ics. Each individual junction is governed by the SNIC 
bifurcation to limit cycle oscillation. We distinguish the 
RCSJ units as excitable when they are in a stable steady 
state for a selected constant current bias less than a crit¬ 
ical value and oscillatory when biased by a higher con¬ 
stant current to cross the SNIC bifurcation point. As a 
result, we find that the presence of a fraction of excitable 
units generates bursting in the whole network although 
the uncoupled oscillatory junctions never show bursting 
dynamics. For a coupling above a threshold, the whole 
network starts synchronous firing with single spiking, and 
for further increase of coupling, periodic bursting appears 
with increasing number of spikes in a single burst and fi¬ 
nally which clearly emerges as a parabolic bursting. Dur¬ 
ing the spiking and bursting above a coupling threshold, 
the whole network splits into two synchronous clusters, 
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one forming a synchronization manifold of the excitable 
units and another of the oscillatory units, however, they 
are phase-locked. We reduce the network model using 
the two synchronization manifolds of the oscillatory and 
excitable units and explain the bursting mechanism and 
furthermore, numerically verify the bursting dynamics of 
the whole network. 


II. SINGLE JUNCTION MODEL 

A single RCSJ model is described by, 

0 + aO + sinO = I. ( 1 ) 

where 0 is phase difference of the junction, 0=v is the 
voltage across the junction, a=[h/2ireIR 2 C J 1 / 2 is the 
damping parameter, h is the Planck’s constant, e is the 
electronic charge and I is a constant bias current. It 
has an equilibrium solution of sin 0=1$ in a cylindrical 
space. The stability of the equilibrium is obtained from 
the /'(#*) = cosO* = (1 — Jq ) 1 / 2 where f'=df/dO at equi¬ 
librium 0 = 0*. For Jo < 1.0, the model has clearly two 
equilibrium points, a node for f r (0 *) < 0 and a saddle 
for f(0 *) > 0. They coalesce at Jo = 1.0 via SNIC bifur¬ 
cation Hi for a choice of a > 1.19. For a < 1.19, a fold 
bifurcation is recorded at Jo = 1.0. In addition there is a 
bistable region for Iq < 1.0 and a < 1.19. We focus here 
on the SNIC regime for Iq > 1.0 and a > 1.19, where 
the stable equilibrium is separated from the oscillatory 
regime by a bifurcation line (Iq = 1.0). 
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FIG. 1. (Color online) Bursting dynamics in a network of 
Josephson junctions. Temporal dynamics shown in the left 
panels and spatio-temporal dynamics in the right panels for 
N=100. Fraction of excitable units in the network, p/N = 0.5. 
I e = 0.5, A = 1.5, p = 0.5,0! — 1.5. Asynchronlous network 
for e = 3.7 (panels in the uppermost row), two synchronous 
clusters with bursting in the second row (e = 5.0), third row 
(e — 8.0), bottom row (e — 9.7). 


III. NETWORK OF JUNCTIONS 


We consider a population of N globally coupled RCSJ 
units in which p number of oscillators are in excitable 
mode (I e . < 1.0), in general, and (N — p) units are self- 
oscillatory (J s . > 1.0). The network consists of two sub¬ 
populations and its dynamics is described by two sets of 
equations, 

N 

6 a + OL ei 0 ei + sinO ei = I ei + — y^($s 7 - — 0 ei ). (2) 

3 = 1 
N 

0 Si + OL Si b Si + sinO Si = I s . + — : ^2(0 e . - 0 S .). (3) 

j =i 



where e* = 1, 2, ...,p and Si = p + 1 ,p + 2,..., N denote 
the excitable and self-oscillatory units respectively. The 
a = 1.5 is chosen identical for all the oscillators to re¬ 
strict our current study in the SNIC regime PQ 3^. The 
bias currents to the excitable and oscillatory units are 
assumed as, J e . = 0.5 and I Si = 1.25 respectively. 

For numerical simulations, we first consider a network 
of size N=100 with equal number of oscillatory and ex¬ 
citable units. Initial conditions are chosen carefully using 
random numbers generated between 0.2 and 0.3. Figure 


FIG. 2. Bifurcation diagram of a Josephson junction unit in 
the network. One oscillatory unit randomly chosen from the 
whole population and shown its bifurcation in the upper panel 
(a), and dynamics of the reduced model at lower panel (b). 
Ie = 0.5, A = 1.5, a = 1.5. 
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FIG. 3. Circle/Circle bursting in the network of Josephson 
junctions. Temporal dynamics in the left panel, inter-spike in¬ 
terval showing parabolic nature in the right panel. Percentage 
of excitable units, p = 0.5, e=20, I e = 0.5, I s = 1.5, a = 1.5. 

[l] reveals a sequence of bursting oscillation in the whole 
network for increasing coupling strength in the upper to 
the lower panels except the uppermost panels. The pan¬ 
els in the uppermost row show no phase-locking for cou¬ 
pling strength e = 3.7. For e > 3.7 in rest of the panels, 
the whole population forms two clusters as seen from the 
time series plot of all the oscillators {0 s i, Si = 50 and 0 e i, 
ei = 50) in each panel. The excitable and the oscillatory 
units form two separate clusters above a threshold cou¬ 
pling (panels in lower three rows), and the two subgroups 
are also seen phase locked. In fact, the first phase-locked 
firing in the whole network starts with single spiking dy¬ 
namics (not shown here) above a coupling threshold and 
then appears the bursting for larger coupling strength 
and adds on one after another spike in each burst (left 
panels in lower three rows). The number of spikes could 
be even larger for further increase of coupling strength as 
shown later, in the text, when we are able to recognize the 
parabolic nature of the bursting. Each of the right pan¬ 
els describes a temporal pattern of all the oscillator nodes 
(si = 50, ei = 50); lower three panels clearly show forma¬ 
tion of two clusters. These are in perfect match with the 
nature of the time series at their immediate left panels. 
This allows a reductionism approach [18j [26] to the large 
network dynamics and restrict them into two synchro¬ 
nization manifolds, 0\ = 6^ = .... = 0 P representing the 
original excitable units and 6 P +1 = Q p +2 = .... = On rep¬ 
resenting the original oscillatory units when we represent 
the network by two oscillators, 

0 e + a e 0 e + sin0 e = I e + e(l — p)(0 s — 0 e ). (4) 

0 S + a 3 0 s + sin0 s = I S + ep{0 e - 0 S ). (5) 

where p/N denotes the fraction of excitable junctions in 
the whole population. 

Figure [2] presents the bifurcation diagram of the dy¬ 
namics of a single oscillatory unit arbitrarily chosen from 
the whole network and its reduced model (4)-(5) as well. 
Maxima of 0 s i of the junction node (say, i = 1) is plotted 
with coupling strength (e) in the upper panel which repre¬ 
sents the original oscillatory units (s^). It shows periodic 
bursting with the number of spikes increasing in a burst 



FIG. 4. (Color online) Circle/Circle bursting of the Josephson 
junctions. Upper panel shows the bursting in a 3D plane of 
0 e , 0 S , Q s . Lower panel shows an enlarged picture of a part 
(in a box) of the upper panel, p — 0.5, e=20, I e = 0.5, I s = 
1.5, a = 1.5. 


one after another with coupling strength. Each period¬ 
adding regime is intercepted by a complex bursting win¬ 
dow. The maxima of 6 S of the reduced model is shown in 
the lower panel and its bifurcation is in agreement with 
the upper panel. The windows of complex dynamics are 
also found matching, which also shows complex bursting 
pattern but here we do not focus on this feature here. 
The reduced model thereby perfectly represents the dy¬ 
namics of the whole network. The excitable units {e/j 
also show similar bifurcation diagram (not shown here) 
and match with the reduced model of the excitable units 
as expected since they are all phase-locked with the os¬ 
cillatory units (s^. 

The nature of bursting is parabolic {circle/circle type) 
mm as shown in Fig. [3j in the sense, that the oscilla¬ 
tion stops via SNIC bifurcation and the trajectory moves 
towards the steady state that becomes unstable via SNIC 
after an elapse of time to start the oscillation once again 
and the process repeats. Left panel shows the time se¬ 
ries of few bursts only; a larger coupling strength is con¬ 
sidered here, when the number of spikes is reasonably 
large. The inter-spike interval in a burst is plotted for 
sucessive spikes in the right panel, which confirmed the 
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parabolic nature of the interspike intervals. The burst¬ 
ing dynamics periodically switches between the oscilla¬ 
tory state and the steady state both via SNIC bifurca¬ 
tion. The excitable units induces a slow dynamics in 
the individual units that controls the firing (oscillatory) 
and steady (resting) states which is further elaborated in 
Fig. B The upper panel demonstrates the circle/circle 
bursting in a 3D plane of 0 e , 0 S , 0 S where the zero plane 
is drawn in cyan/gray color. A part of the trajectory 
(in a box) goes below the zero plane which is enlarged 
in the lower panel. The location of the saddle with its 
eigen-directions and the node of the uncoupled excitable 
junction are denoted by open and soild circles respec¬ 
tively. 

We further explain the bursting mechanism here as 
controlled by the slow dynamics. The trajectory of the 
bursting, after a few fast spiking, moves towards the sad¬ 
dle (plane in cyan/gray line) guided by the stable man¬ 
ifold of the excitable unit. The trajectory obviously be¬ 
comes slow when it approaches the saddle point, and after 
coming sufficiently close to it, moves away by the influ¬ 
ence of the unstable egienvector. We make an approxi¬ 
mation here that the excitable and the oscillatory units 
behave like isolated units since all the oscillators are now 
phase locked, when the error dynamics is negligbly small; 
the coupling term in the system is almost negligible. The 
ensemble of exctibale units and oscillatory units are re¬ 
duced to two junctions. The excitable unit maintains its 
isolated dynamics with a node and a saddle that influ¬ 
ences the slow-fast bursting of the oscillatory junction. 
As a result, the network of junctions show a bursting 
dynamics typical of the Type I neuron. In fact, we ob¬ 
serve the circle/circle bursting simply by coupling two 
oscillators, one excitable and another oscillatory. 


IV. CONCLUSION 

In summary, we investigated a mixed population of os¬ 
cillatory and excitable Josephson junctions under all-to- 
all global coupling when we observed circle/circle burst¬ 


ing in a broad parameter range of the junction and the 
coupling strength. We produced numerical evidence of 
the phenomenon using a network of N=100 oscillators 
and taking two equal populations of oscillatory and ex- 
citbale junctions. The whole network splits into two 
clusters for our chosen range of coupling strength that 
helps reduce the system into a two-oscillator model. Re¬ 
sults of the reduced model were found perfectly match¬ 
ing with the numerical results of the whole network. We 
found that the number of spikes increases with coupling 
strength which we supported with a bifurcation diagram 
of the whole network and its reduced model. The burst¬ 
ing dynamics had been a dominant feature of the mixed 
population such that it existed for different percentage 
of excitable units although we have only deatiled the 
case of fifty-fifty populations of oscillatory and excitable 
units. The bursting is typically circle/circle type in the 
selected parameter regime of the superconducting device 
where the dynamics is governed by the SNIC bifurcation. 
We have numerically checked that the bursting also ap¬ 
peared in a mixed population of another SNIC model, 
the Morris-Lecar system [25] . The bursting dynamics, 
therefore, seemed to be a generic feature of a mixed pop¬ 
ulation of such dynamical units and of even a larger size. 
A mixed population of oscillatory and excitable units, 
in general, was investigated earlier by others [18, 119], 
however, they focussed on the global oscillation, partial 
oscillation and oscillation death regimes. In the glob¬ 
ally synchronized (1:1 and higher phase locking) oscilla¬ 
tion regime, the type of dynamics was mentioned m as 
could be of complex type, however, it was not focussed. 
Our results, in that sense, provided additional informa¬ 
tion about the network dynamics of a mixed population 
in the global oscillation regime (1:1 phase locking) and, 
particularly, showed clear evidence of bursting in a SNIC 
model, the superconducting junction. Other bifurcation 
regimes of the junction such as the fold bifurcation and 
the bistable region are also interesting which we plan to 
explore, in the future, especially using other complex net¬ 
work topologies. 
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